!> \file GFS_suite_interstitial_3.F90
!!  Contains code to setup convectively-transported tracers, calculate critical relative humidity, and save cloud number concentrations

  module GFS_suite_interstitial_3

  contains

!> \section arg_table_GFS_suite_interstitial_3_run Argument Table
!! \htmlinclude GFS_suite_interstitial_3_run.html
!!
    subroutine GFS_suite_interstitial_3_run (otsptflag,                 &
               im, levs, nn, cscnv,imfshalcnv, imfdeepcnv,              &
               imfshalcnv_samf, imfdeepcnv_samf,progsigma,              &
               first_time_step, restart,                                &
               satmedmf, trans_trac, do_shoc, ltaerosol, ntrac, ntcw,   &
               ntiw, ntclamt, ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc,    &
               xlon, xlat, gt0, gq0, sigmain,sigmaout,qmicro,           &
               imp_physics, imp_physics_mg,                             &
               imp_physics_zhao_carr, imp_physics_zhao_carr_pdf,        &
               imp_physics_gfdl, imp_physics_thompson, dtidx, ntlnc,    &
               imp_physics_wsm6, imp_physics_fer_hires, prsi, ntinc,    &
               imp_physics_nssl,                                        &
               prsl, prslk, rhcbot,rhcpbl, rhctop, rhcmax, islmsk,      &
               work1, work2, kpbl, kinver, ras, me, save_lnc, save_inc, &
               ldiag3d, qdiag3d, index_of_process_conv_trans,           &
               clw, rhc, save_qc, save_qi, save_tcp, errmsg, errflg)

      use machine, only: kind_phys

      implicit none

      ! interface variables
      logical, intent(in)     :: otsptflag(:)!  on/off switch for tracer transport (size ntrac)
      integer,              intent(in   )                   :: im, levs, nn, ntrac, ntcw, ntiw, ntclamt, ntrw, ntsw,&
        ntrnc, ntsnc, ntgl, ntgnc, imp_physics, imp_physics_mg, imp_physics_zhao_carr, imp_physics_zhao_carr_pdf,   &
        imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6,imp_physics_fer_hires,  &
        imp_physics_nssl, me, index_of_process_conv_trans
      integer,              intent(in   ), dimension(:)     :: islmsk, kpbl, kinver
      logical,              intent(in   )                   :: cscnv, satmedmf, trans_trac, do_shoc, ltaerosol, ras, progsigma
      logical,              intent(in   )                   :: first_time_step, restart
      integer,              intent(in   )                   :: imfshalcnv, imfdeepcnv, imfshalcnv_samf,imfdeepcnv_samf 
      integer,                                          intent(in) :: ntinc, ntlnc
      logical,                                          intent(in) :: ldiag3d, qdiag3d
      integer,              dimension(:,:),             intent(in) :: dtidx
      real,                 dimension(:,:),            intent(out) :: save_lnc, save_inc

      real(kind=kind_phys), intent(in   )                   :: rhcbot, rhcmax, rhcpbl, rhctop
      real(kind=kind_phys), intent(in   ), dimension(:)     :: work1, work2
      real(kind=kind_phys), intent(in   ), dimension(:,:)   :: prsl, prslk
      real(kind=kind_phys), intent(in   ), dimension(:,:)   :: prsi
      real(kind=kind_phys), intent(in   ), dimension(:)     :: xlon, xlat
      real(kind=kind_phys), intent(in   ), dimension(:,:)   :: gt0
      real(kind=kind_phys), intent(in   ), dimension(:,:,:) :: gq0

      real(kind=kind_phys), intent(inout   ), dimension(:,:)   :: sigmain
      real(kind=kind_phys), intent(inout   ), dimension(:,:)   :: sigmaout,qmicro
      real(kind=kind_phys), intent(inout), dimension(:,:)   :: rhc, save_qc
      ! save_qi is not allocated for Zhao-Carr MP
      real(kind=kind_phys), intent(inout), dimension(:,:)   :: save_qi
      real(kind=kind_phys), intent(inout), dimension(:,:)   :: save_tcp
      real(kind=kind_phys), intent(inout), dimension(:,:,:) :: clw

      character(len=*),     intent(  out)                   :: errmsg
      integer,              intent(  out)                   :: errflg

      ! local variables
      integer :: i,k,n,tracers,kk
      real(kind=kind_phys) :: tem, tem1, tem2
      real(kind=kind_phys), dimension(im) :: tx1, tx2, tx3, tx4

      !real(kind=kind_phys),parameter :: slope_mg = 0.02, slope_upmg = 0.04,  &
      !                   turnrhcrit = 0.900, turnrhcrit_upper = 0.150
      ! in the following inverse of slope_mg and slope_upmg are specified
      real(kind=kind_phys), parameter :: zero = 0.0_kind_phys, one = 1.0_kind_phys
      real(kind=kind_phys), parameter :: slope_mg   = 50.0_kind_phys,   &
                                         slope_upmg = 25.0_kind_phys

      ! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

      ! In case of using prognostic updraf area fraction, initialize area fraction here
      ! since progsigma_calc is called from both deep and shallow schemes.
      if(((imfshalcnv == imfshalcnv_samf) .or. (imfdeepcnv == imfdeepcnv_samf)) &
           .and. progsigma)then
         if(first_time_step .and. .not. restart)then
            do k=1,levs
               do i=1,im
                  sigmain(i,k)=0.0
                  sigmaout(i,k)=0.0
                  qmicro(i,k)=0.0
               enddo
            enddo
         endif
         do k=1,levs
            do i=1,im
               sigmaout(i,k)=0.0
            enddo
         enddo
      endif


      if (cscnv .or. satmedmf .or. trans_trac .or. ras) then
        tracers = 2
        do n=2,ntrac
!          if ( n /= ntcw  .and. n /= ntiw  .and. n /= ntclamt .and. &
!               n /= ntrw  .and. n /= ntsw  .and. n /= ntrnc   .and. &
!               n /= ntsnc .and. n /= ntgl  .and. n /= ntgnc) then
            IF ( otsptflag(n) ) THEN
            tracers = tracers + 1
            do k=1,levs
              do i=1,im
                clw(i,k,tracers) = gq0(i,k,n)
              enddo
            enddo
          endif
        enddo
      endif ! end if_ras or cfscnv or samf

      if (ntcw > 0) then
        if (imp_physics == imp_physics_mg .and. rhcpbl < 0.5_kind_phys) then ! compute rhc for GMAO macro physics cloud pdf
          do i=1,im
            tx1(i) = one / prsi(i,1)
            tx2(i) = one - rhcmax*work1(i)-rhcbot*work2(i)

            kk     = min(kinver(i), max(2,kpbl(i)))
            tx3(i) = prsi(i,kk)*tx1(i)
            tx4(i) = rhcpbl - rhctop*abs(cos(xlat(i)))
          enddo
          do k = 1, levs
            do i = 1, im
              tem  = prsl(i,k) * tx1(i)
              tem1 = min(max((tem-tx3(i))*slope_mg, -20.0_kind_phys), 20.0_kind_phys)
              ! Using rhcpbl and rhctop from the namelist instead of 0.3 and 0.2
              ! and rhcbot represents pbl top critical relative humidity
              tem2 = min(max((tx4(i)-tem)*slope_upmg, -20.0_kind_phys), 20.0_kind_phys) ! Anning
              if (islmsk(i) > 0) then
                tem1 = one / (one+exp(tem1+tem1))
              else
                tem1 = 2.0_kind_phys / (one+exp(tem1+tem1))
              endif
              tem2 = one / (one+exp(tem2))

              rhc(i,k) = min(rhcmax, max(0.7_kind_phys, one-tx2(i)*tem1*tem2))
            enddo
          enddo
        else
          do k=1,levs
            do i=1,im
              kk = max(10,kpbl(i))
#ifdef SINGLE_PREC
              if (k < kk) then
                tem    = rhcbot - (rhcbot-rhcpbl) * (one-prslk(i,k)) / max(one-prslk(i,kk),1e-7)
              else
                tem    = rhcpbl - (rhcpbl-rhctop) * (prslk(i,kk)-prslk(i,k)) / max(prslk(i,kk),1e-7)
              endif
#else
              if (k < kk) then
                tem    = rhcbot - (rhcbot-rhcpbl) * (one-prslk(i,k)) / (one-prslk(i,kk))
              else
                tem    = rhcpbl - (rhcpbl-rhctop) * (prslk(i,kk)-prslk(i,k)) / prslk(i,kk)
              endif
#endif
              tem      = rhcmax * work1(i) + tem * work2(i)
              rhc(i,k) = max(zero, min(one,tem))
            enddo
          enddo
        endif
      else
        rhc(:,:) = 1.0
      endif

      if (imp_physics == imp_physics_zhao_carr .or. imp_physics == imp_physics_zhao_carr_pdf) then   ! zhao-carr microphysics
        !GF* move to GFS_MP_generic_pre (from gscond/precpd)
        ! do i=1,im
        !   psautco_l(i) = Model%psautco(1)*work1(i) + Model%psautco(2)*work2(i)
        !   prautco_l(i) = Model%prautco(1)*work1(i) + Model%prautco(2)*work2(i)
        ! enddo
        !*GF
        do k=1,levs
          do i=1,im
            clw(i,k,1) = gq0(i,k,ntcw)
          enddo
        enddo
      elseif (imp_physics == imp_physics_gfdl) then
        clw(1:im,:,1) = gq0(1:im,:,ntcw)
      elseif (imp_physics == imp_physics_thompson) then
        do k=1,levs
          do i=1,im
            clw(i,k,1)    = gq0(i,k,ntiw)                    ! ice
            clw(i,k,2)    = gq0(i,k,ntcw)                    ! water
            save_tcp(i,k) = gt0(i,k)
          enddo
        enddo
        if(ltaerosol) then
          save_qi(:,:) = clw(:,:,1)
          save_qc(:,:) = clw(:,:,2)
        else
          save_qi(:,:) = clw(:,:,1)
        endif
      else if (imp_physics == imp_physics_nssl ) then
        do k=1,levs
          do i=1,im
            clw(i,k,1) = gq0(i,k,ntiw)                    ! cloud ice
            clw(i,k,2) = gq0(i,k,ntcw)                    ! cloud droplets
          enddo
        enddo
          save_qi(:,:) = clw(:,:,1)
          save_qc(:,:) = clw(:,:,2)
      elseif (imp_physics == imp_physics_wsm6 .or. imp_physics == imp_physics_mg .or. imp_physics == imp_physics_fer_hires) then
        do k=1,levs
          do i=1,im
            clw(i,k,1) = gq0(i,k,ntiw)                    ! ice
            clw(i,k,2) = gq0(i,k,ntcw)                    ! water
          enddo
        enddo
      endif

      if(imp_physics == imp_physics_thompson .and. ldiag3d .and. qdiag3d) then
         if(dtidx(100+ntlnc,index_of_process_conv_trans)>0) then
            save_lnc = gq0(:,:,ntlnc)
         endif
         if(dtidx(100+ntinc,index_of_process_conv_trans)>0) then
            save_inc = gq0(:,:,ntinc)
         endif
      endif

    end subroutine GFS_suite_interstitial_3_run

  end module GFS_suite_interstitial_3
